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I propose a large class of stochastic Markov processes associated with probability distributions analogous to 
that of lattice gauge theory with dynamical fermions. The construction incorporates the idea of approximate 
spectral split of the determinant through local loop action, and the idea of treating the infrared part of the split 
through explicit diagonalizations. I suggest that exact algorithms of practical relevance might be based on the 
Markov processes so constructed. 



One of the difficulties we face in lattice gauge 
theory is how to efficiently simulate the gauge in- 
variant distribution 



let P 2 occ 



(U, U') = min[l, T^jjf] be the Metropo- 
lis acceptance probability with respect to P 2 . 
Then the transition matrix 



P{U) oc e~ SB{u) det M(U) = P 1 (U)P 2 (U) (1) T(U, U') = T X (U, f/')P 2 acc (f/, U% U ^ U' (2) 



with bounded bosonic action Sb{U), and the 
fermionic kernel M(U), such that det M(U) > 0. 

The dynamic Monte Carlo approach, which is 
usually adopted, is based on identifying a suitable 
ergodic Markov process, defined by the Markov 
matrix T(U,U'), which has required distribution 
as a fixed point, i.e. / dUP(U)T(U, U') = P{U'). 
While finding examples of valid Markov processes 
is by no means difficult, finding a suitable one 
proved to be a formidable task. 

In the discussion that follows, I will construct 
Markov processes built on rather different ingre- 
dients than those of established methods. They 
incorporate several interesting ideas that ap- 
peared in recent years, and the corresponding ex- 
act algorithms have a chance of being reasonably 
practical, while perhaps being more accomoda- 
tive of complicated actions of Ginsparg- Wilson 
type. The main influence in the construction I 
will present is the method used in the inexact 
Truncated Determinant Algorithm (TDA) 0. 

Probability distribution that splits into two 
unnormalized products such as ([l) can in prin- 
ciple be simulated using the following state- 
ment. Let Ti(U,U') be the Markov matrix sat- 
isfying detailed balance with respect to Pi, i.e. 
Pi(U)Tx(U,U')dU = PiiU^T^U'^dU', and 
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satisfies detailed balance with respect to the dis- 
tribution P oc P1P2. Consequently, T represents 
a valid Markov process assuming T\ is ergodic. 

While there is a lot of freedom in choosing the 
process Ti, the algorithm based on T very much 
exemplifies the difficulties associated with simu- 
lating distribution (Q). For one, the calculation 
of acceptance probability requires the calculation 
of determinant ratio, which is hard. Also, if T± 
represents the probabilities for the sweep of suit- 
able local changes, for example, then the determi- 
nant is expected to fluctuate strongly, giving rise 
to severe problems with equilibration and, hence, 
useless algorithm. Nevertheless, this is a natu- 
ral starting point and the aim is to improve upon 
these bad things. 

As a first point to notice, there is an additional 
significant freedom in the above procedure. In 
particular, it is possible to split the distribution 
(Q) in two parts in infinitely many ways by multi- 
plying and dividing with arbitrary bounded pos- 
itive function of the gauge fields, i.e. 



P 1 (U) = 



-Sb-ASe 



P 2 (U) = e ASB detM (3) 



Labeling the splits by ASb, we thus have a fam- 
ily of Markov processes T[ASb] constructed as in 
(||), each with the correct fixed point. 

The reason why this can be useful is implicitly 
contained in the well established fact, that if the 
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fermions are not very light, the effects of the de- 
terminant can be absorbed in the simple shift of 
the gauge coupling . This suggests that one can 
follow the fluctuations of the determinant in such 
a case by approximating it with a simple loop ac- 
tion. Using such an action for the split ASb as 
described above can therefore presumably reduce 
the fluctuations of the determinant part substan- 
tially even if the fermions are not heavy. Analo- 
gous observations in more or less related contexts 
were made for example in Refs. |]-[|]. 

Supposing this works, we are still left with 
an awkward algorithm, because the evaluation of 
Metropolis acceptance probability is costly. One 
way to proceed is to consider the stochastic linear 
Kennedy-Kuti acceptance probability instead ||] , 
which can be cheap for the price of introduc- 
ing some extra noise. Indeed, we can general- 
ize the transition matrices (||) to stochastic ones 
and use the following statement. Let Ti(U,U') 
be the Markov matrix satisfying detailed balance 
with respect to Pi. Let further R a {U,U') be 
a stochastic estimator, depending on the set of 
stochastic variables a, so that < R a (U,U') > a = 
P2(U')/P 2 (U), and let 

pace = (\+ + X-R a (U,U'), if U>U'; U) 
2 ' a \A_ + X+R a (U,U'), if U<U', K > 

where A + , A_ are constants and some ordering of 
the gauge fields is assumed. If we define 

f a (U, U') = T^U, U')P^ a c (U, U'), U^U' (5) 

then T(U, U') =< f a (U, U') > a satisfies detailed 
balance with respect to P oc P\P 2 . Consequently, 
the stochastic transition matrices (|^) can serve to 
generate a valid Markov sequence assuming that 
T\ is ergodic and the individual estimates P 2 C £ 
can be interpreted as probabilities. 

With this framework in mind, we are now 
dealing with large class of stochastic matrices 
f[AS B ,P^ a c ] that are assigned to P(U) and, in 
addition to the choice of split action, are also la- 
beled by the choice of stochastic probability esti- 
mator. The challenge now is to construct an esti- 
mator that is (a) reasonably cheap, (b) does not 
require too much noise so that the fluctuations 
already reduced by the split will not come back, 



(c) introduces negligible amount of probabilistic 
violations in the acceptance step. 

That this can be done in the useful novel way 
is suggested by the qualitative observation that 
the small eigenvalues of the lattice Dirac opera- 
tor contribute substantially to the long distance 
behaviour of the effective fermionic action while 
large eigenvalues are more relevant at short dis- 
tances . This obviously makes very good phys- 
ical sense, but it should be said that to make 
more quantitative statements appears to be non- 
trivial. Nevertheless, it suggests that after the 
fluctuations of the full determinant are reduced 
by splitting the distribution with some ultralocal 
loop action, what is left in the determinant part 
is dominated by small eigenvalues. We should 
thus construct an estimator R a {U,U') entering 
(|J) that assumes bulk contribution from small 
eigenvalues, thus being able to take advantage of 
the approximate spectral split. 

Such an estimator can be constructed as fol- 
lows. By definition, we need to estimate the quan- 
tity R(U,U') = P 2 {U')/P 2 {U). If Ai denote the 
eigenvalues of M(U), this can be written as 



as b {u>) X [X' 2 ...X' N _ 
e AS B (u) \ X X 2 ...X 



XiX 2 . . . X N , 



N 



(6) 



where Xi 



A, 



e [AS B (U')-AS B (U)]/N ^ and the 

eigenvalues are assumed to be ordered by increas- 
ing magnitude. The monomial can then be esti- 
mated by 



Ra 



X\ + (l\ — (x 2 — 1) - 

aia 2 {x 3 - 1) + 



(7) 



a\a 2 



+ aia 2 . 



X\X 2 . . . X N -i 

■ ajv-i {xn 



1) 



a±a 2 . 



, aN-i 



where < ai < 1 and ai are independent random 
variables distributed according to 



1, with probability a*; 
0, with probability 1 — ai 



(8) 



Obviously, one has < R a > a = Xix 2 . . . x n as de- 
sired, and it is useful to emphasize the following 
points: 

(A) It is trivial to generate di by using the ran- 
dom number generator. To get i? Q , one can first 
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assign to it the value x\ and generate a\. If a\ is 
zero, this is the whole estimate because a\ mul- 
tiplies all the remaining terms. If a± is one, the 
second term is added and 0,2 is generated. Again, 
if a,2 is zero, the estimate is completed and if it 
is one, the third term is added and &3 generated 
etc. The calculation stops as soon as the first a,i is 
zero, and only the corresponding number of small- 
est eigenvalues is needed. The required number 
of eigenvalues can be determined beforehand if 
desired, and the average number over many esti- 
mates can be tuned by changing {c^}. 

(B) One can group the contributions of several 
eigenvalues into a single variable X{ , and the num- 
ber of eigenvalues so grouped does not have to be 
the same for all Xi. For example, if the split action 
ASb is determined by fitting the ultraviolet part 
of the truncated determinant Q , then it is proba- 
bly best to group into x\ the smallest eigenvalues 
almost up to the truncation number, and group 
the remaining eigenvalues differently so that the 
stochastic part of R a does not fluctuate a lot. 

(C) Strictly speaking, \ do not have to be the 
eigenvalues of M. For example, in case of two 
flavours of Wilson fermions, we can instead use 
Xi = S*5i, where Si are the eigenvalues of one 
flavour operator. It is in fact desirable that the 
corresponding Xi be real non-negative and so, if 
the individual eigenvalues are not, we assume that 
we can group them together so that the resulting 
product is, or use additional properties such as in 
the example above. This is possible in situations 
of practical interest. 

(D) The amount of noise introduced by the es- 
timator (Q) is typically much less than with the 
traditional e TrlogM estimators, and can be tuned 
by changing the values {ai}. This comes at the 
price of calculating the lowest eigenvalues. 

For obvious reasons, it is natural to call the al- 
gorithms based on f[AS B , P2.T] of (§) with the 
estimator of type (^) the Stochastic Split Deter- 
minant Algorithms (SSDA). They simultaneously 
use the complementary representations of effec- 
tive fermionic action in terms of gauge loops with 
increasing length (converging rapidly at small dis- 
tances), and in terms of Dirac eigenvalues with 
increasing magnitude (converging rapidly at large 
distances) . If the two representations overlap suf- 



ficiently with small number of terms, then a sim- 
ple split action ASb can be found and T\, 
easily adjusted, so that the number of proba- 
bilistic exceptions is negligible, and efficient ex- 
act SSDA results. The TDA work (|§ suggests 
that this might be the case for reasonably large 
lattices in QCD4. The work on quantitative as- 
pects of these statements is in progress and will 
be reported elsewhere. 

On sufficiently large lattices, algorithms like 
TDA or SSDA in their current form will eventu- 
ally have inferior efficiency to that of HMC, which 
has more favourable volume scaling. However, 
the crucial advantage of TDA and SSDA is that 
they essentially treat fermions in the eigenspace 
of operator M. As such, they are in principle ap- 
plicable to functions f(M) as well. For example, 
simulations of two or arbitrary flavours of stag- 
gered fermions is straightforward here. 

Finally, it should be emphasized that in the 
current context, the underlying process T\ can 
always be chosen so that there will essentially 
be no exceptions at all. In bad cases this can 
rapidly increase the work per independent con- 
figuration, but exactness can always be achieved. 
However, the underlying philosophy appears to 
be sufficiently solid to believe that the framework 
is large enough for some practically relevant al- 
gorithms to be found here. 

I thank Hank Thacker for many pleasant dis- 
cussions on the topics presented here. 
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